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ABSTRACT 

The  requirement  that  the  propellants  used  in  launch  vehicle 
systems  must  also  provide  the  thermal  energy  to  be  converted 
to  kinetic  energy  in  the  rocket  nozzle  has  plagued  rocket 
designers  since  the  dawn  of  the  space  age.  Beamed  propulsion 
systems,  however,  avoid  this  constraint  by  placing  the  energy 
source  on  the  ground  and  transmitting  the  energy  to  the 
spacecraft  via  microwaves.  This  work  computationally  models 
three  different  heat  exchanger  channel  designs  for  use  in  a 
beam  propulsion  spacecraft.  It  was  found  that  despite  the  very 
small  diameter  of  the  microchannels,  each  design  produced 
extreme  temperature  gradients  across  the  channel  cross  section. 

INTRODUCTION 

Currently  the  choice  of  propellants  used  in  launch  vehicles 
is  limited  by  the  requirement  that  the  propellants  also  must 
provide  the  thermal  energy  that  is  to  be  converted  into  kinetic 
energy  in  the  nozzle.  Because  of  this  requirement,  modern  bi¬ 
propellant  launch  systems  must  carry  both  a  fuel  and  an 
oxidizer  such  as  in  the  case  of  the  space  shuttle,  which  carries 
both  hydrogen  and  oxygen.  In  addition  to  decreasing  the 
specific  impulse  of  the  rocket  below  that  of  a  system  using 
hydrogen  as  its  only  propellant,  one  must  also  carry  separate 
tanks,  pumps,  and  plumbing  if  an  oxidizer  is  to  be  carried  in 
addition  to  the  fuel.  For  the  past  five  decades  bi-propellant 
rocket  propulsion  has  been  the  singular  method  available  to 
launch  large  payloads  into  space.  Because  of  this,  bi-propellant 


rocket  designs  have  already  been  improved  and  refined 
significantly.  Therefore,  a  new  method  of  rocket  propulsion  is 
required  if  any  significant  gains  in  rocket  performance  are  to  be 
achieved  [1]. 

The  concept  of  beamed  energy  propulsion  eliminates  the 
need  for  the  propellant  to  provide  its  own  thermal  energy. 
Rather  than  increasing  the  temperature  of  a  gas  through 
chemical  reactions  between  a  fuel  and  an  oxidizer,  the  energy 
source  is  kept  on  the  ground  and  transmitted  to  the  launch 
vehicle  through  the  use  of  microwaves  or  lasers.  A  wide  variety 
of  different  beamed  energy  propulsion  concepts  exist  and  the 
mm-wave  energy  sources  required  are  technologically  feasible. 
However,  in  order  for  this  concept  to  work,  a  heat  exchanger 
must  be  designed  that  can  efficiently  transmit  the  heat  given  by 
the  energy  into  the  hydrogen  propellant  [2].  This  work 
computationally  models  a  proposed  heat  exchanger  tube  for  this 
purpose. 

Much  work  has  been  done  on  the  design  and  effectiveness 
of  microchannel  heat  exchangers  for  the  purposes  of  cooling 
electronics  [3].  However,  these  heat  exchangers  work  at  far 
lower  temperatures  than  those  that  are  of  interest  here.  The 
ultimate  goal  of  this  work  is  to  design  and  characterize  a  heat 
exchanger  that  efficiently  absorbs  the  incoming  mm-waves  and 
couples  the  thermal  energy  to  the  hydrogen  flow.  Temperatures 
exceeding  2000  K  are  expected.  This  paper  presents  a 
preliminary  feasibility  study  for  a  high  temperature  hydrogen 
heat  exchanger.  Only  the  coupling  of  the  thermal  energy  to  the 
flow  is  investigated.  Microchannels  of  three  different  shapes 
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are  evaluated  and  the  resulting  temperatures  and  velocities  from 
the  numerical  model  are  presented. 

NOMENCLATURE 

Tmax  =  maximum  temperature 
Vmax  =  maximum  velocity 
Tave  =  mass  flow  weighted  average  temperature 
AP  =  pressure  drop 
Acs  =  cross-sectional  area 
Arad  =  area 
Dh  =  hydraulic  diameter 
q"  =  heat  flux  per  unit  area 
a  =  Stefan-Boltzmann  constant 

MODEL 

Using  the  COMSOL  Multiphysics  4.3  package,  three- 
dimensional  models  of  all  three  different  heat  exchanger  tube 
designs  were  created.  A  cross-sectional  view  of  the  three  tube 
shapes  is  shown  below  Figure  1.  All  three  shapes  were  sized  to 
have  a  hydraulic  diameter  of  1  mm.  Each  tube  is  .  1  m  long  and 
made  entirely  of  the  working  fluid.  As  it  was  assumed  that  the 
material  that  formed  the  tube  walls  had  an  infinite  thermal 
conductivity  the  solid  tube  itself  was  not  modeled.  The  only 
effect  caused  by  this  assumption  was  that  the  heat  flux  was 
applied  directly  to  the  hydrogen  rather  than  to  the  surrounding 
solid  tube. 

In  order  to  analyze  the  specified  situation,  the  heat  transfer 
in  fluids  and  low  Reynolds  number  k-s  turbulent  flow  physics 
modules  were  coupled  together.  This  method  was  found  to  be 
more  stable  than  the  use  of  COMSOL’s  conjugate  heat  transfer 
physics  module.  The  two  physics  modules  used  were  coupled 
simply  by  setting  the  heat  transfer  physics  to  use  the  fluid 
pressure  and  velocity  of  the  turbulent  flow  module  and  setting 
the  turbulent  flow  module  to  use  the  temperature  given  by  the 
heat  transfer  module. 

Each  model  consisted  of  a  single  tube  of  hydrogen  that  was 
governed  by  a  set  of  six  basic  boundary  conditions.  Three  of 
these  constrained  the  heat  transfer  felt  by  the  hydrogen  while 
the  other  three  defined  the  fluid  flow.  Thus  all  of  the  tube’s 
boundaries  were  constrained  by  two  boundary  conditions;  one 
for  the  heat  transfer  physics  module  and  one  for  the  turbulent 
flow  physics  module. 

For  the  inflow  portion  of  the  tube  both  a  constant  flow 
velocity  and  a  constant  temperature  boundary  condition  were 
applied.  The  magnitude  of  this  velocity  was  dictated  by  a  mass 
flow  rate  requirement  of  4.66E-5  kg/s.  Set  to  match  the  initial 
temperature  of  hydrogen  as  it  exits  from  its  storage  tank,  the 
constant  temperature  condition  prescribed  a  temperature  of  300 
K  at  the  inlet. 

For  the  walls  of  the  tube  a  no-slip  boundary  condition  was 
applied  in  conjunction  with  a  heat  flux  condition.  The  value  of 
the  heat  flux  was  a  two  part  calculation.  First,  the  incoming 
heat  flux  was  calculated  so  that  for  each  shape  a  total  flux  of 
1.33  kW  would  be  achieved.  Next,  from  this  calculated  input 
value  the  amount  of  heat  lost  due  to  outgoing  radiation  was 
subtracted.  Here,  based  on  prospective  tube  shell  materials 


such  as  carbon  fiber  reinforced  carbon,  an  emissivity  of  .9  was 
chosen.  The  resultant  temperature  dependent  equation  is  shown 
below  as  Equation  1.  In  this  equation,  a  is  the  Stefan- 
Boltzmann  constant  and  T  is  the  temperature  of  the  wall  at  the 
point  where  this  equation  is  applied.  This  equation  was  used  for 
the  cylindrical  model.  In  the  other  two  models,  the  same  basic 
equation  was  used  but  the  incoming  heat  flux  term  was 
modified  so  that  design  would  receive  a  total  of  energy  input  of 
1.3  kW  before  radiative  effects  are  included. 

q"  =  4.23352L6  -  (.9  *  a  *  ( T 4  -  3004))  (1) 

Lastly,  to  the  exit  of  the  tube  both  heat  outflow  and 
pressure-no  viscous  stress  boundary  conditions  were  applied. 
To  the  latter  of  these  conditions  a  gage  pressure  of  0  kPa  was 
set. 

Finally,  to  save  time  on  computation,  symmetry  boundary 
conditions  were  used  for  the  square  and  triangular  cross- 
sectional  tubes.  For  example,  the  square  tube  was  cut  along  its 
lines  of  symmetry  so  that  only  a  quarter  of  the  square  remained. 
Symmetry  boundary  conditions  were  then  applied  to  the  sides 
along  the  lines  of  symmetry.  This  method  enabled  a 
significantly  greater  mesh  resolution  for  the  same  amount  of 
computational  time.  Figure  1  shows  the  triangular  and  square 
tubes  before  and  after  symmetry  conditions  were  used. 


Figure  1  -  Tube  cross  sections  after  applying  symmetry 
conditions.  Here,  the  solid  lines  show  the  final  shape,  the 
dotted  lines  show  the  removed  portions  and  the  red  dashed 
lines  show  the  lines  of  symmetry  used. 

After  the  boundary  conditions  had  been  applied  to  the 
system,  two  adjustments  were  made  to  the  default  physics 
settings.  First,  the  reference  absolute  pressure  in  both  modules 
was  set  to  be  1724  kPa,  the  pressure  given  by  the  preliminary 
design.  Second,  the  heat  transfer  module  was  set  to  use  the 
ideal  gas  law  to  calculate  the  density  of  the  hydrogen. 

While  COMSOL’s  material  library  can  provide  some  of  the 
needed  material  properties  of  hydrogen,  certain  temperature 
dependent  values,  such  as  thermal  conductivity  and  dynamic 
viscosity,  are  only  given  by  COMSOL  for  temperatures  up  to 
1300  K.  Because  of  this,  COMSOL’s  built  in  hydrogen  material 
data  was  modified  to  include  material  properties  for  the  higher 
temperatures  experienced  by  the  heat  exchanger  designs 
evaluated  here  [4]. 

For  the  purposes  of  these  models  it  was  assumed  that  the 
hydrogen  remained  entirely  in  its  diatomic  form  and  did  not 
experience  any  dissociation  due  to  the  high  temperatures  it  was 
heated  to.  While  the  validity  of  this  assumption  begins  to  break 
down  around  the  maximum  temperatures  experienced  by  the 
hydrogen  in  these  models,  due  to  the  limited  amounts  of  the 
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hydrogen  experiencing  these  maximum  temperatures  and 
previous  research  showing  that  the  amount  of  atomic  -  rather 
than  diatomic  -  hydrogen  should  be  far  less  than  7%  [5],  the 
amount  of  error  resulting  from  this  assumption  has  been 
evaluated  to  be  negligible. 

Rather  than  using  COMSOL’s  default  mesh,  a  custom 
mesh  was  created  for  each  tube.  This  was  done  in  three  steps. 
First,  the  inlet  face  of  the  modeled  tube  was  meshed  with  a 
triangular  mesh.  The  resolution  of  the  mesh  was  increased  at 
the  outer  walls  of  each  tube  in  order  to  provide  greater  detail  of 
the  boundary  layer.  The  second  step  was  to  sweep  the  mesh 
from  the  inlet  face  to  the  outlet  face.  Finally,  the  boundary 
layer  mesh  feature  was  also  added  to  the  tube  to  further 
increase  resolution  at  the  tube’s  outer  walls.  This  custom 
meshing  resulted  in  meshes  of  236800,  244400,  and  88000 
elements  for  the  square,  triangular,  and  circular  tubes, 
respectively.  The  lower  total  element  number  of  the  cylindrical 
model  is  due  to  the  fact  that  the  axisymmetric  model  was  able 
to  provide  much  more  detail  for  far  less  computational  time. 
Thus  the  axisymmetric  model  was  used  for  analysis  of  the 
cylindrical  heat  exchanger  and  the  cylindrical  three  dimensional 
model  was  used  for  comparison  only. 

The  default  iterative  solver  chosen  by  COMSOL  proved  to 
be  the  most  efficient  in  obtaining  a  solution  for  the  three- 
dimensional  models.  Instead  of  attempting  to  solve  for  each  of 
the  dependent  variables  at  the  same  time,  the  solution  process 
was  broken  up  into  three  different  segregated  steps.  In  each  step 
the  iterative  solver  was  used  in  conjunction  with  a  multigrid 
algorithm  rather  than  an  Incomplete  LU  algorithm  in  order  to 
reduce  computational  time. 

VALIDATION 

In  addition  to  the  three-dimensional  tubes,  a  two- 
dimensional  axisymmetric  model  of  the  cylindrical  tube  was 
created.  The  creation  of  this  model  served  several  purposes. 
First,  by  utilizing  the  axial  symmetry  of  the  cylindrical  tube  in 
the  axisymmetric  model  the  tube  could  be  treated  as  a  two- 
dimensional  -  rather  than  a  three-dimensional  -  problem.  This 
allowed  for  an  extreme  increase  in  mesh  resolution  for  the  same 
computational  time  as  the  three-dimensional  model.  Second,  by 
evaluating  the  cylindrical  tube  using  both  the  mesh  from  the 
three-dimensional  models  and  the  axisymmetric  model  it  was 
possible  to  prove  that  all  of  the  solutions  obtained  by  COMSOL 
for  this  particular  analysis  were  mesh  independent.  Finally, 
computing  both  models  allowed  for  validation  of  the  method  of 
using  symmetry  boundary  conditions  as  was  applied  to  the 
square  and  triangular  tubes.  Significant  data  from  the 
axisymmetric  model  is  given  below  in  Table  2. 


Table  1  -  Significant  data  from  the  axisymmetric  model 


Tmax  [K] 
Vmax  [m/s] 
Tave[K] 
AP  [kPa] 
Acs  [mm2] 
AJmm2] 
Dh  [mml 


Axisymmetric 

2482.6 

353.4 

1588 

36.3 

.785 

314.16 

1 


RESULTS 

Despite  the  geometric  similarity  of  the  models,  the  results 
from  the  numerical  analysis  showed  some  significant 
differences.  Significant  values  from  each  model  are  presented 
below  in  Table  2. 


Table  2  -  Si; 

gnificant  data  from  each  of  the  4  mode 

Square 

Triangle 

Cylinder 

Tmax  [K] 

2590.7 

2580 

3470 

Vmax  [m/s] 

276 

213 

356 

Tave  [K] 

1541 

1438 

1602 

AP  [kPa] 

24.2 

15.6 

43.5 

Acs  [mm2] 

1 

1.3 

.785 

AJmm2] 

400 

519.6 

314.16 

Dh  [mm] 

1 

1 

1 

Figure  2,  below,  shows  the  velocity  profile  at  the  exit  of  the 
square  tube.  Unexpectedly,  the  maximum  velocity  did  not  occur 
in  the  middle  of  the  tube  (that  is,  the  top  right  of  this 
symmetrical  model)  but  rather  was  located  in  between  the 
tube’s  center  and  outer  wall.  This  is  explained  by  the  fact  that 
this  is  a  temperature  driven  flow  rather  than  a  pressure  driven 
flow.  In  a  fully-developed  pressure  driven  flow  the  farther  one 
proceeds  from  the  tube  wall  across  the  cross  section  of  the  tube 
the  higher  the  flow’s  velocity  gets.  However,  in  each  of  the 
models  presented  here  the  driving  force  behind  the  flow  is 
temperature  rather  than  pressure.  Thus,  while  the  frictional 
effects  experienced  by  pressure  driven  flows  are  still  present 
here,  their  effect  on  the  velocity  profile  is  partially  countered  by 
the  nature  of  the  temperature  driven  flow  found  here. 
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Figure  2  -  Velocity  profile  at  exit  plane  of  square  model. 
Note  location  of  maximum  velocity. 

As  the  temperature  of  the  hydrogen  increases,  the  density 
decreases.  As  the  density  of  the  hydrogen  decreases,  the  flow 
velocity  will  increase  assuming  a  constant  cross-sectional  area, 
as  shown  by  the  conservation  of  mass  in  Equation  2. 


positioned  between  the  tube’s  center  and  its  outer  wall. 
However,  unlike  the  square  tube  where  the  maximum  velocity 
occurs  at  four  points  around  the  tube  center,  the  maximum  fluid 
velocity  in  the  triangular  tube  occurs  at  three  distinct  points 
within  the  tube. 


Fitil  Vnlacriy  [m/s| 
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Figure  4  -  Velocity  of  hydrogen  at  the  exit  of  the  triangular 
tube. 


P\A1v1  —  P2A2V2  (2) 

Thus  if  frictional  effects  are  ignored,  the  hottest  part  of  the 
flow  will  also  be  the  fastest.  As  Figure  3  shows,  the  hottest  part 
of  the  flow  is  at  the  tube  walls.  The  combination  of  the  wall 
friction  and  the  temperature  driven  nature  of  the  flow  produces 
a  velocity  profile  in  which  the  point  of  highest  velocity  is 
neither  at  the  center  nor  at  the  walls  of  the  tube,  but  rather  lies 
somewhere  in  between  these  two  extremes. 


fcatrt  l  enp erasure  LK| 
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Figure  3  -  Temperature  at  exit  of  square  model. 


Figure  5  provides  the  temperature  of  the  hydrogen  at  the 
exit  of  the  triangular  tube.  As  expected  the  tips  of  the  triangle 
experienced  the  maximum  temperature.  While  it  is  expected 
that  if  the  cross-sectional  area  of  each  model  had  been  held 
constant  the  maximum  temperature  of  the  triangular  model 
would  have  been  higher  than  that  of  the  square  model,  due  to 
the  fact  the  hydraulic  diameter  of  each  model  was  chosen  to  be 
held  constant  rather  than  the  area  the  square  model  possess  the 
highest  maximum  temperature. 

fcxit  l  enp erasure  IN 
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Figure  5  -  Temperature  at  the  exit  of  the  triangular  model. 


The  results  from  the  triangular  model  were  similar  to  the 
results  from  the  square  model.  As  Figure  4  shows,  the 
maximum  fluid  velocity  within  the  hydrogen  is  once  again 


Figures  6  and  7,  respectively,  present  the  exit  velocity  and 
temperature  of  the  three-dimensional  cylindrical  model. 
However,  as  the  axisymmetric  model  offers  the  same 
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Exit  velocity  [nvsj 


information  as  the  three-dimensional  model  but  with  a  much 
higher  accuracy  due  to  the  axisymmetric  model’s  far  greater 
resolution,  all  further  discussion  of  results  of  the  cylindrical 
tube  shape  will  be  restricted  to  the  axisymmetric  model. 

FhM,  ViJuctiy  [m/s| 
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Figure  6  -  Fluid  velocity  at  the  exit  plane  of  the  three- 
dimensional  cylindrical  model 


•30  -20  -10 


10  20  30  40  50  xlO*  ▼  -4. 1397x10’*' 


Figure  8  -  Fluid  velocity  at  the  exit  boundary  of  the 
axisymmetric  model 
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Figure  7  -  Fluid  temperature  at  the  exit  plane  of  the  three- 
dimensional  cylindrical  model 

Figures  8  and  9  depict  the  exit  velocity  and  temperature  as 
calculated  by  the  axisymmetric  model.  It  is  clearly  seen  that 
these  figures  closely  parallel  Figures  6  and  7.  This  similarity 
both  confirms  the  validity  of  methods  of  symmetry  used  here 
and  also  verifies  the  mesh  independent  nature  of  the  solutions 
obtained. 
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Figure  9  -  Fluid  temperature  at  the  exit  boundary  of  the 
axisymmetric  model 

Figures  10  and  11,  respectively,  give  the  cross-sectional 
temperature  and  velocity  of  the  fluid  as  a  function  of  position 
along  the  tube  length  as  obtained  by  the  axisymmetric  model. 
In  Figure  10,  the  flat  line  at  0  m  shows  the  set  300  K 
temperature  at  the  inlet  boundary.  The  temperature  of  the 
system  increases  with  position  along  the  tube.  The  wall 
temperature  (width  =  50E-5  m)  increases  more  rapidly  than 
flow  core  (width  =  0  m)  temperatures.  Half  way  down  the  tube, 
the  wall  temperature  is  2300  K  but  the  tube  center  is  still  at  300 
K.  At  the  tube  exit,  the  wall  temperature  reaches  2500  K  but  the 
flow  center  is  at  550  K.  The  flow  velocities  and  wall  effects 
produce  this  large  gradient.  The  wall  temperature  increases  are 
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most  dramatic  near  the  tube  entrance,  increasing  1500  K  in  the 
first  0.1  m.  Beyond  the  tube  entrance,  changes  in  wall 
temperature  are  significantly  smaller  due  to  the  decreased 
gradient. 

Although  a  constant  inlet  velocity  boundary  condition  was 
used  for  the  inlet  boundary  conditions  on  each  of  the  three- 
dimensional  models  due  to  the  relative  instability  of  the 
constant  mass  flow  boundary  condition,  due  to  the  comparative 
simplicity  of  the  axisymmetric  model  the  more  accurate 
constant  mass  flow  boundary  condition  was  used.  Because  of 
this,  the  velocity  across  the  entrance  -  specified  as  the  0  m  line 
in  Figure  11  -  is  not  constant  but  rather  shows  a  more  accurate 
variance  between  the  center  of  the  tube  (width  =  0  m)  and  the 
tube  wall  (width  =  50E-5  m).  This  plot  clearly  shows  that 
everywhere  along  the  tube  length,  except  the  tube  entrance,  the 
maximum  velocity  does  not  occur  at  the  tube  center,  but 
somewhere  in  between  the  tube  center  and  wall.  The  location  of 
the  maximum  moves  towards  the  tube  center  as  the  tube  exit  is 
neared.  This  can  be  attributed  to  the  increased  penetration  of 
the  heating  into  the  flow  with  tube  length  as  seen  in  Figure  10. 

Figures  12  and  13  show  the  temperature  and  velocity 
profiles  along  the  length  of  the  tube.  In  each  case,  the  line 
marked  R  =  0  m  gives  the  profile  at  the  center  of  the  tube  and 
the  line  marked  R  =  .0005  m  gives  the  profile  at  the  tube  wall. 
As  expected,  Figure  12  shows  that  the  fluid  temperature 
increased  fastest  and  became  the  highest  along  the  walls  of  the 
tube.  The  wall  temperature  exceeds  2000  K  for  over  80%  of  the 
tube  length.  In  Figure  13,  it  can  be  seen  that  the  velocity  trends 
change  beyond  a  tube  length  of  0.01  m.  Near  the  entrance,  the 
flow  acts  like  a  typical  pressure  driven  flow  with  the  slowest 
flow  at  the  wall  and  the  velocity  increasing  towards  the  center. 
Beyond  the  entrance,  the  flow  velocity  quickly  reaches  a 
maximum  at  a  radius  of  2E-4  m.  Clearly,  the  effect  of  heating  is 
felt  by  the  flow  almost  immediately. 
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Figure  10  -  Hydrogen  temperature  across  cross-section  of 
axisymmetric  model 
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Figure  11  -  Fluid  velocity  across  cross-section  of 

axisymmetric  model 
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Figure  12  -  Hydrogen  temperature  along  the  length  of  the 
tube 
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Figure  13  -  Fluid  velocity  along  the  length  of  the  tube 


CONCLUSION 

While  the  models  analyzed  in  this  paper  may  be  simple, 
the  numerical  study  provided  several  interesting  results.  First, 
temperature  gradients  of  almost  2000  K  were  seen  across  the 
radius  of  each  of  the  three  tubes.  These  extreme  gradients  were 
not  expected  due  the  tube’s  microscale  radii.  Secondly,  the 
velocity  profiles  produced  across  the  cross  section  of  each 
shape  were  not  the  expected  paraboloids.  Rather  than  showing 
the  previously  expected  exponential  rise  in  velocity  from  its 
minimum  at  the  tube  wall  to  its  maximum  at  the  tube  center, 
each  model’s  velocity  profile  placed  the  maximum  fluid 
velocity  at  a  location  in  between  the  tube  center  and  wall. 
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